Damage-cluster distributions and size effect on strength in compressive failure 
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We investigate compressive failure of heterogeneous materials on the basis of a continuous pro- 
gressive damage model. The model explicitely accounts for tensile and shear local damage and 
reproduces the main features of compressive failure of brittle materials like rocks or ice. We show 
that the size distribution of damage-clusters, as well as the evolution of an order parameter, the size 
of the largest damage-cluster, argue for a critical interpretation of fracture. The compressive failure 
strength follows a normal distribution with a very small size effect on the mean strength, in good 
agreement with experiments. 



Understanding how materials break is a fundamen- 
tal problem that has both theoretical and practical rele- 
vance. The topic has received considerable renewed at- 
tention during the last few decades because of the limi- 
tations of the classical Griffith theory for heterogeneous 
media [l|43|. The practical applications are numerous, 
from material and structural design, to the important 
problem of size effects on strength (e.g. [J]). The two 
key components that make material failure challenging 
to understand arc long-range interactions and material 
disorder. 

Traditionally, Wcibull and Gumbel distributions asso- 
ciated with the weakest-link approach have been widely 
used to describe the strength of brittle materials. These 
distributions naturally arise from extreme- value statistics 
of initial defect distributions based on the assumptions 
that J5( (i) defects do not interact with one another, (ii) 
failure of the whole system is dictated by the activation of 
the largest flaw (the weakest- link hypothesis) , and finally 
(iii) the material strength can be related to the critical 
defect size. These assumptions are reasonable for mate- 
rials with relatively weak disorder loaded under tension, 
but do not hold for heterogeneous materials with broad 
distribution of initial disorder or for loading conditions 
stabilizing crack propagation, such as compression. In 
these cases there is experimental evidence that a consid- 
erable amount of damage can be sustained before failure 
p. Under these conditions, failure is the culmination 
of a complex process involving the nucleation, propa- 
gation, interaction and coalescence of many microcracks 
[7J. Stress states observed under various natural condi- 
tions, ranging from soil and rock mechanics to earthquake 
physics, suggest the importance of compressive failure. 

A cornerstone for the understanding of breakdown of 
disordered media has been lattice models of fracture in 
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which networks with prescribed bond failure thresholds 
are subject to increasing external loads [2|- In these mod- 
els, failure is described on a qualitative level as the inter- 
play between disorder and elasticity. When strong dis- 
order is considered, these models suggest that fracture 
strength does indeed not follow a Weibull or a Gumbel 
distribution but a log- normal distribution [8J, |9( . Similar 
strength distributions have been obtained for different 
model types (fuses and springs), in 2D and 3D, suggest- 
ing that it is a general feature of failure in heterogeneous 
materials with broad disorder [2|. 

On a more fundamental point of view, the evolution 
of the distribution of crack-cluster sizes in the vicinity 
of failure has implications on the interpretation of rup- 
ture as a phase transition. For a second-order or criti- 
cal phase transition, local quantities such as the size of 
crack-clusters are expected to show scaling and a diverg- 
ing correlation length, such as in the percolation problem 
[l|. In the limit of infinitely strong disorder, fracture can 
be mapped onto the percolation problem [lfj , suggesting 
a critical transition. This interpretation remains how- 
ever controversial in the case of non-infinite disorder, as 
lattice models show an abrupt localization of damage at 
failure, without a diverging correlation length, arguing 
instead for a first-order transition |9|, LLD, [LSI ■ 

In this letter, we revisit these problems for compressive 
failure of an heterogeneous material with variable range 
of disorder, on the basis of a continuous progressive dam- 
age model [Tj, [Tj] ■ This model is more realistic than lat- 
tice models or scalar damage models [15| , as it explicitely 
accounts for the tensorial nature of stresses and strains. 
The principal features of compressive failure of brittle 
materials like low-porosity rocks or ice are captured by 
the model: the macroscopic stress-strain response or the 
progressive localization of damage onto an inclined fault 
at failure [g, |lj| . A detailed comparison of the model 
with experimental results of compressive failure of rock 
samples has been performed using acoustic emission and 



damage avalanche statistics [l7j. Recently, an analysis 
of this damage localization, either tracked from damage 
avalanches or from the evolution of continuous strain-rate 
fields, showed a divergence of the associated correlation 
length towards the failure, i.e. argued for the critical 
point hypothesis [lj|. Here we report that the size dis- 
tribution of damage clusters as well as the evolution of 
an order parameter, the size of the largest damage clus- 
ter, also argue for a critical interpretation of compressive 
fracture with specific scaling laws in the pre- as well as 
post-failure phases. We also show that the compressive 
failure strength has a normal distribution and character- 
ize the associated size effect. 

The model, described in more details elsewhere [1J, 
[jjjl, considers a continuous 2D clastic material (Hooke's 
law) under plane stress, with progressive local damage. 
Damage is represented by a reduction of the isotropic 
elastic modulus Yi of the element i, Yi(n + 1) = Yi(n)do 
with do = 0.9, each time the stress state on that element 
exceeds a given threshold. This elastic softening simu- 
lates an increase in crack density at the element scale 
[18[ as supported by experiments [lj|. In high poros- 
ity materials, modifications of elastic properties can also 
result from local compaction (e.g. [20|), a problem not 
considered here. The stress field is recalculated each time 
a damage event occurs by solving the equation of static 
equilibrium using a finite element scheme. As the result 
of elastic interactions, the stress redistribution following 
a damage event can set off an avalanche of damage, which 
stops when the damage threshold is no longer fulfilled by 
any element. 

The Coulomb criterion, r = jjlun + C, of wide ap- 
plicability for brittle materials under compressive stress 
states [21] . defines the damage threshold, r and ctjv 
are respectively the shear and normal stress on the el- 
ement (sign convention positive in compression), fi is 
an internal friction coefficient identical for all elements, 
whereas quenched disorder is introduced through the co- 
hesion C randomly drawn from a uniform distribution. 
We use /i = 0.7, a common value for most geomaterials 
|2l| . This envelope is completed by a truncation in ten- 
sion in the Mohr's plane, i.e. the clement is damaged if 
cjv = &Ntensiie = — 2xl0~ 3 x>o- The simulations, which 
start with undamaged material (Yi = Yq = const.), are 
performed on rectangular meshes of randomly oriented 
triangular elements. A uniaxial compression loading is 
applied by increasing the vertical displacement of the up- 
per boundary of the system (i.e. strain-driven loading), 
whereas the lower boundary is fixed and left and right 
boundaries can deform freely. The loading increment is 
extrapolated to damage the weakest element, ensuring 
infinitely slow driving compared to stress redistribution 
time. Two series of simulations, with meshes of linear size 
L varying from 8 to 128 elements, were performed with 
different ranges of disorder: 0.5 x 10~ 3 F < C < 10~ 3 F 
which we refer to as the HI disorder, 0.2 x 10 _3 Yo < C < 
10 _3 Yn refered to as H2. The range of disorder used in 
this study is thus narrower than in lattice models where 
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FIG. 1. (a) Example of the evolution of the macroscopic stress 
versus strain, as well as the damage avalanche size, i.e. num- 
ber of damage events per iteration (L = 64, H2 heterogene- 
ity). Red arrows indicate respectively (i) peak load and (ii) 
the point where the largest damage cluster spans the system. 
Maps of damage clusters at these two points are respectively 
plotted as (b) and (c) with the largest cluster in red. 



the failure threshold distribution usually extends down 
to (e.g. J8|). The number of independent simulations 
performed with each system size is 10 4 up to L = 16, 
5 x 10 3 for L = 32, 10 3 for L = 64 and 10 2 for L = 128. 

In the early stages of deformation, damage is scat- 
tered homogeneously and the macroscopic stress-strain 
response remains essentially linear (Fig. [Ua)). As defor- 
mation proceeds, macroscopic softening occurs up to a 
maximum stress, the strength or peak load, followed by 
one or a few macroscopic stress drops. Then, the macro- 
scopic stress remains approximately constant, i.e. the 
behaviour mimics plasticity. 

First, we identify damage clusters whose size A is de- 
fined as the total surface area of adjacent elements (i.e. 
sharing two nodes) that have been damaged at least one 
time from the beginning of the simulation. At peak load, 
the largest damage cluster does not yet span the system 
(Fig. H{b)). Then, during the post-peak phase, damage 
events are localized in the vincinity of one or a few large 
damage clusters which eventually evolve into a spanning 
cluster, connecting two opposite boundaries of the sys- 
tem (Fig. HJc)). Using similar simulations, we previously 



reported [lj] that the spatial correlation length associ- 
ated with damage events or strain localization reaches 
the size of the system at peak load. This means that the 
divergence of the dynamical correlation length associated 
with damage avalanches and strain localization precedes 
the geometrical evolution of the fault. In other words, 
due to the long-range correlations in the stress field that 
emerge in the vicinity of failure, damage clusters do not 
need to be interconnected into a spanning cluster for the 
global failure to occur. Strain-driven compressive failure 
experiments of rocks have reported similar observations: 
the failure plane is not fully formed at peak load and only 
appears in the post-peak phase with acoustic emissions 
strongly localized along the main fault during this phase 
|6j. This agreement with experiments indicates that the 
model is relevant to describe damage evolution in the 
post-peak phase. 

In what follows the distance to the peak load, identi- 
fied as the critical point, is tracked in terms of a control 
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where e m is the macro- 



scopic strain, reaching e mp at peak load. This definition 
yields A = 1 at the first damage event, decreasing down 
to A = at the peak load and increasing again after peak 
load. 

We analyse the distribution of damage cluster sizes A: 
before peak load it follows a power law with an expo- 
nential cut-off at large sizes that increases in the vicinity 
of peak load (Fig. |2J), P(A) ~ A~ a &xp(-A/A*). For 
finite-size systems, the evolution of A* is controled by 
the system-size itself (finite size effect) as well as the dis- 
tance to peak load. We make the simplifying hypothesis 
that, for the largest system (L = 128), the cut-off size 
is only a function of the control parameter, A* ~ A r . 
This allows us to estimate the value of a and 7 through a 
data collapse (Fig. [21 inset). For the HI disorder we find 
a = 3.6±0.1 and 7 = 1.3±0.2, while for H2, a = 2.6±0.1 
and 7 = 1.6 ±0.2. 

The nature of the distribution identified argues for a 
critical interpretation of failure, since we show here that 
a local quantity, the damage cluster size, shows scaling 
in the vicinity of failure. Moreover, this evolution is re- 
lated to the divergence of the correlation length and the 
critical exponents can be related geometrically: since A is 
the surface area of damage clusters, it can be expected to 
scale with the correlation length £ and the fractal dimen- 
sion of damage clusters D as A* ~ £ D . In Ref . [1J] we re- 
ported from an analysis of strain-rate fields £ ~ A A , with 
A = 1.0 ± 0.1, D = 1.15 ± 0.05 for HI, and D = 1.4 ± 0.1 
for H2. The geometrical relationship A* ~ A XD yields 
7 = XD for large system sizes. When accounting for the 
uncertainty on the exponent values, the later relationship 
is captured for the largest system size (L = 128), with 
a correlation length exponent independent of disorder. 
The analysis can also be refined to account for the fi- 
nite size effect on A [22[ . These results contrast with the 
log-normal distribution of crack-cluster sizes observed for 
lattice models of fracture |23| . 

After peak load, the distribution of cluster sizes shows 
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FIG. 2. (a) Probability density function (PDF) of the size of 
damage clusters A for bins of the control parameter A before 
peak load, and associated data collapse (inset), (b) PDF of 
A after the peak load (L = 128, H2 disorder). The color code 
corresponds to the value of A as given by the legend, red is 
closest to peak load. 



a power law with similar exponent, as well as an ad- 
ditional peak at large sizes that progressively grows as 
the distance to failure increases (Fig. [5]). This suggests 
that the physics governing the evolution of damage is 
different during this phase and favors the growth of the 
largcsts clusters. As the changes in the distribution are 
concentrated in its tail during this phase, we analyze the 
growth of the largest damage cluster, which size is de- 
noted n. By analogy with the percolation theory, we 
consider n as the order parameter of the phase transi- 
tion, defined here for the entire deformation history. The 
general evolution of n is continuous with an inflexion 
point at failure (Fig. 131(a)). The degree of initial hetero- 
geneity shows a large influence on the behavior around 
peak load, a narrower heterogeneity leading to a more 
abrupt transition. At peak load, II scales with the sys- 
tem size U(L, A ->■ 0) = U p (L) = L s , with S = 0.3 ± 0.1 
for HI and 5 = 0.8 ± 0.1 for H2. This means that with 
a broader disorder the largest cluster is in average closer 
to span the system at peak load. 

Next, we analyze the growth of the largest cluster in 
the post-peak phase. The net growth of n scales with 
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FIG. 3. (a) Growth of the largest damage cluster size IT, 
A = corresponds to peak load, left and right of this point are 
respectively the pre- and post-peak phases, (b) Data collapse 
of the largest cluster size II after the peak load (H2 disorder). 
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FIG. 4. Normal probability plot of the strength distribution. 
Albeit minute deviation, the plot shows a collapsed straight 
line, as expected for a normal distribution (L = 128, H2). The 
figure also demonstrates the power law scaling of the standard 
deviation std(<Tf) ~ L~ v . 



the system size as well as with the control parameter 
following: 



(n.-n. p )~L 5 'A . 



(i) 



The value of 6' is estimated from a data collapse, yielding 
6' = 1.3 ± 0.1 (HI) and 5' = 1.9 ± 0.1 (H2) (Fig. 1 (b)). 
In the vicinity of peak load, the resulting collapsed curve 
shows a power law yielding /3 — 0.4 (HI) and /3 = 0.7 
(H2). 

We established that the growth of clusters during this 
post-peak phase is equally contributed by the coalescence 



of existing clusters as well as by the extension of the total 
damaged area, since the two mechanisms are described 
by a scaling law equivalent to Eq. [H with similar expo- 
nents |22| ■ The two mechanisms are tightly linked during 
this phase, as the expansion/branching of clusters bridges 
gaps with other clusters and induces coalescence. Close 
to peak load, we observed that Eq. Q] also applies to 
the growth of other large damage clusters, i.e. the 2nd, 
3rd, 4th largests clusters for example, until they even- 
tually coalesce or stop growing as distance to peak load 
increases. 

On Fig. [3] (b), the deviation from the power law for 
A > 0.1 corresponds to the point where the largest clus- 
ter becomes spanning. Beyond this point the geometrical 
growth of the largest cluster appears progressively inhib- 
ited and the power law scaling does not hold, i.e. further 
deformation is accomodated only by additional damage 
and shear along the inclined, mature fault, in full agree- 
ment with observations |6j. 

Finally, we analyse the compressive failure strength, 
which we define as the maximal macroscopic stress cr/ 
reached during a simulation. We verified that the 
strength distribution cannot be represented by a Wcibull 



nor by a Gumbcl distribution [22j, meaning that the 
weakest-link approach is irrelevant to compressive fail- 
ure. Instead, the strength distribution can be described 
by a normal distribution (Fig. 01 . Fig. 2] also shows 
that the standard deviation of the strength scales as 
std(<jf) ~ L~ v , where a data collapse is obtained for 
r\ = 0.4 for HI disorder and r/ = 0.65 for H2. The mean 
strength shows a very slow decrease with increasing sys- 
tem size that can be represented as (<jf) ~ L~ e where 
6 = 0.02. This is in excellent agreement with compres- 
sive failure experiments on rock and ice that show no 
significant sample size effect on the mean strength [24J . 
To conclude, we showed that in compressive failure, 
the peak load can be considered as a critical point re- 
garding the evolution of damage clusters. Scaling laws, 
specific to the pre- and post-peak phases, were shown to 
describe the evolution of the size distribution of damage 
clusters. This expresses the difference in the physics that 
govern the growth of damage clusters in these two phases. 
Compressive failure thus appears as a complex cumu- 
lative process involving long-range correlations, interac- 
tions and coalescence of microcracks. In such conditions, 
the hypotheses of the weakest-link approach, describing 
failure as an abrupt fist-order transition, are violated 
and the strength distributions predicted by extreme- 
value statistics do not apply at all. In lattice models 
of fracture, strength distributions are neither captured 
by Weibull nor by Gumbel distributions [3, Q. An in- 
teresting point is that the ranges of initial disorder that 
we have considered are not as broadly distributed as in 
lattice models. This suggests that for loading modes sta- 
bilizing damage propagation, such as compression, even 
under narrowly distributed initial disorder, the weakest- 
link approach docs not hold and failure can be interpreted 
as a critical phase transition. 



All computations were performed at SCCI-CIMENT Grenoble. 
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